source("utils.R")
source("gamm4_utils.R")
load_all_pkgs() 
theme_set(theme_bw())
zmargin <- theme(panel.spacing=grid::unit(0,"lines"))
zmarginx <- theme(panel.spacing.x=grid::unit(0,"lines"))
library('pander')  ## for tables
sapply(pkgList,function(x) format(packageVersion(x)))
##          lme4         gamm4         bbmle   broom.mixed          brms 
## "1.1.21.9002"       "0.2.5"    "1.0.23.1"  "0.2.4.9000"      "2.11.1" 
##    data.table       lattice     gridExtra       ggplot2       viridis 
##      "1.12.8"     "0.20.38"         "2.3"       "3.2.1"       "0.5.1" 
##        plotly       cowplot      ggstance          plyr         dplyr 
##       "4.9.1"       "1.0.0"       "0.3.3"       "1.8.5"       "0.8.4" 
##         tidyr        tibble         remef        r2glmm        raster 
##       "1.0.2"       "2.1.3"  "1.0.6.9000"       "0.1.2"      "3.0.12" 
##         rgdal        fields       plotrix            sp 
##       "1.4.8"        "10.3"       "3.7.7"       "1.3.2"
comb_out <- function(p,fn,...) {
    print(p)
    htmlwidgets::saveWidget(ggplotly(p,...),fn)
}
best_models <- readRDS("bestmodels_gamm4.rds")
ecoreg <- readRDS("ecoreg.rds")
ecoreg <- ecoreg

Plots of the original variables:

taxon_names <- c("mbirds","mamph","mmamm")
respvars <- paste0(taxon_names,"_log")
predvars <- c("NPP_log_sc","Feat_log_sc","NPP_cv_sc","Feat_cv_sc")
## plot raw values of diversity for a given taxon/predictor
plotfun(model=NULL, respvar="mbirds_log",
        predvar="NPP_log_sc",
        backtrans=TRUE,log="xy")

## generate all combinations
orig_plots <- list()
for (i in seq_along(taxon_names)) {
    orig_plots[[taxon_names[i]]] <- list()
    for (pp in predvars) {
      orig_plots[[taxon_names[i]]][[pp]] <-
        plotfun(model=NULL,respvar=respvars[i],
                xvar=pp,data=ecoreg,
                backtrans=TRUE,log="xy")+
                theme(legend.position="none")
    }
}
plot_grid(plotlist=orig_plots[["mbirds"]],align="hv",nrow=2)

Plots of the transformed variables (log and scaled)

## names are the names of the variable in the data set;
##  the character strings are the names to put on the axis
## this time using log-scaled etc.
pred_names <- c(NPP_mean="NPP g*C/m2*year",
                Feat_mean="Feat (% of NPP)",
                NPP_cv_inter="CV of NPP",
                Feat_cv_inter = "CV of Feat")
sc_taxon_names <- paste(c("Birds","Amphibians","Mammals"),
                      "N sp/km2")
names(sc_taxon_names) <- paste0(c("mbirds","mamph","mmamm"),"_log")
sc_pred_names <- c(NPP_log_sc="log NPP g*C/m2*year",
                   Feat_log_sc="log Feat (% of NPP)",
                   NPP_cv_sc="scaled CV of NPP",
                   Feat_cv_sc = "scaled CV of Feat")

do_plots_1 <- function(taxon="mbirds",pred="NPP_mean",
                       s_names=taxon_names,
                       p_names=pred_names) {
    ggplot(ecoreg,aes_string(x=pred,y=taxon,colour="biome")) +
        geom_point() +
        theme(legend.position="none") +
        labs(y=s_names[taxon], x=p_names[pred])
}

## do all combinations
sc_plots <- list()
for (tt in names(sc_taxon_names)) {
    sc_plots[[tt]] <- list()
    for (pp in names(sc_pred_names)) {
        sc_plots[[tt]][[pp]] <- do_plots_1(tt,pp,
                                           p_names=sc_pred_names,
                                           s_names=sc_taxon_names)
    }
}

## for a given taxon, draw unscaled plots in row 1 and scaled plots in row 2
both_plots <- function(taxon) {
    plot_grid(plotlist=c(orig_plots[[taxon]],
                         sc_plots[[paste0(taxon,"_log")]]),
              nrow=2,
              align="hv")
}

FIXME: y-axis names

Plots of residuals

I generated several “lists” of plots (using plotfun) with the lapply function following Ben’s code. In each case, I changed the xvar and auxvar to plot the 4 top most variables for each taxa (aside from NPP)

taxon <- "mbirds_log"
plotfun(best_models[[taxon]],
        respvar=taxon,
        xvar='Feat_log_sc',auxvar=NULL,data=ecoreg,backtrans=TRUE,log="xy")

## partial residuals plots
remef_plot <- function(taxon="mbirds_log",predvar="NPP_log",
                       auxvar=NULL,title=NULL) {
    m <- best_models[[taxon]]
    if (is.null(title)) {
        title <- if (is.null(auxvar)) predvar else {
             paste(predvar,auxvar,sep=":")                                     
                                              }
    }
    pp <- (plotfun(m,xvar=predvar,respvar="partial_res",
                   auxvar=auxvar,data=ecoreg,
                   backtrans=TRUE,log="xy") 
      + theme(legend.position="none")
      + ggtitle(title)
    )
    return(pp)
}

## rp1: NPP_log
## rp2: NPP_log:Feat_cv_sc
## rp3: NPP_log:Feat_log
## rp4: Feat_log:NPP_cv_sc
## rp5: Feat_log
## rp6: NPP_cv_sc
## rp7: Feat_cv_sc
rem_predvars <- list("NPP_log_sc",
                     "NPP_log_sc",
                     "NPP_log_sc",
                     "Feat_log_sc",
                     "Feat_log_sc",
                     "NPP_cv_sc",
                     "Feat_cv_sc")
rem_auxvars <- list(NULL,"Feat_cv_sc","Feat_log_sc",
                    "NPP_cv_sc",NULL,NULL,NULL)

## construct all combinations of partial residuals for each taxon
remef_plots <- list()
for (tt in names(sc_taxon_names)) { ## for each taxon
    remef_plots[[tt]] <- list()
    for (pp in seq_along(rem_predvars)) { ## for each predictor variable
      ## cat(tt,pp,"\n")
      ## cat(tt," ",pp,"\n")
      ## generate and save the partial residuals plot
      remef_plots[[tt]][[pp]] <- remef_plot(tt,predvar=rem_predvars[[pp]],
                                            auxvar=rem_auxvars[[pp]])
      ## print(remef_plots[[tt]][[pp]])
    }
}

Graphical outputs

Graphic example to see the color of each biome in the plots below

ggplot(ecoreg,aes(NPP_mean,mbirds,colour=biome)) + geom_point() + labs(x = "NPP (g*C/m2*year)", y = "Birds (N sp/km2)")

Amphibians

Original + transformed values (xy plot)

both_plots("mamph")

Top 4 effects

do.call(grid.arrange,
        c(remef_plots[["mamph_log"]][c(2,3,4,5)],
          ## NPP_log:Feat_csv, NPP_log:Feat_log, Feat_log:NPP_cv_sv, Feat_log
          list(ncol=2)))

Birds

Original + transformed values

both_plots("mbirds")

Top 4 effects

do.call(grid.arrange,
        c(remef_plots[["mbirds_log"]][c(5,2,6,7)],
               list(ncol=2)))

Mammals

Original + transformed values

both_plots("mmamm")

Top 4 effects

do.call(grid.arrange,
        c(remef_plots[["mmamm_log"]][c(7,5,2,4)],
          list(ncol=2)))

questions

  • why do fitted model lines look bad? Is this just some kind of problem with the intercept? Under what circumstances do we expect this to work? Scale (geometric mean problem), or incorrect model center, or … ?